Self-consistent Langmuir waves in resonantly driven thermal plasmas 
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The longitudinal dynamics of a resonantly driven Langmuir wave are analyzed in the limit that the 
growth of the electrostatic wave is slow compared to the bounce frequency. Using simple physical ar- 
guments, the nonlinear distribution function is shown to be nearly gaussian in the canonical particle 
action, with a slowly evolving mean and fixed variance. Self-consistency with the electrostatic po- 
tential provide the basic properties of the nonlinear distribution function including a frequency shift 
that agrees well with driven, electrostatic particle simulations. This extends earlier work on nonlin- 
ear Langmuir waves by Morales and O'Neil [G.J. Morales and T.M. O'Neil, Phys. Rev. Lett. 28, 417 
(1972)], and could form the basis of a reduced kinetic treatment of Raman backscatter in a plasma. 

PACS numbers: 52.35.Fp, 52.35. Mw, 52.35.Sb 



43 

Oh! 



CZ3 

a 

11; 
o ■ 

CZ2 . 
>v 

43 



> ■ 

00 ■ 
00 1 
(N . 

o ■ 
r» : 
o . 

o : 

>>', 

43 : 

Or 



X 



I. INTRODUCTION 

Much of the progress in understanding Langmuir waves 
has been from the linear viewpoint, obtained by as- 
suming that the perturbation of the plasma from its 
(Maxwellian) equilibrium is "sufficiently small" such that 
second order terms in the perturbation may be neglected. 
Under these conditions, one can derive the normal modes 
of the distribution function (the singular Case- Van Kam- 
pen modes [H, 0) or, alternatively, the Landau damped 
"modes" of the electric field [3, 4]. A basic result of these 
linear analyses is that smooth, electrostatic perturbations 
tend to zero through the process of Landau damping. 

That such damping is not universal was first pointed 
out by Bernstein, Green, and Kruskal (BGK) [5(, who in- 
cluded the particles trapped in the electrostatic wave to 
formulate nonlinear distribution functions that give rise 
to time-independent electrostatic disturbances. Explicit 
constructions of sinusoidal, small-amplitude BGK waves 
were later derived by Holloway and Doming Q , in which 
they showed that even arbitrarily small amplitude waves 
can exist without being Landau damped. BGK distribu- 
tions are generally functions of the conserved particle en- 
ergy H = ^m e v 2 — e$(z) whose charge perturbation self- 
consistently generates the electrostatic potential &(z) via 
Poisson's equation. Thus, BGK distributions are static 
solutions to the ID Vlasov-Poisson system 



dt 



f{v,z\t) 



df 
dt 



df <9$ df 







(la) 



0_ 2 

dz 



dz dz dv 
4ire I dv f(v, z;t) — 4irerii(z), (lb) 



where / is the electron distribution function and we con- 
sider the ions to have a time-independent density rii(z). 

In this paper, we introduce and characterize nonlin- 
ear Langmuir wave solutions to the Vlasov-Poisson sys- 
tem ([T]) that are naturally occurring BGK-like waves. 
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These waves (and the distribution functions that gen- 
erate them) have particular relevance to laser-plasma 
physics, in that they dynamically arise as kinetic, non- 
linear Langmuir waves in systems that are weakly driven 
on or near resonance. To obtain these solutions, we use 
the canonical action-angle coordinates, finding that the 
plasma is well-described by a simplified distribution func- 
tion that is gaussian in the canonical action. In this way, 
we obtain near-equilibrium solutions that approximate 
the fully time-dependent distribution function when the 
resonant forcing is small. While these notions may be 
reminiscent of adiabatic theory, we do not invoke adi- 
abatic invariance, since we imagine that the plasma is 
resonantly driven. Our calculation is more in the spirit 
of an averaged theory, in that the dynamical dependence 
on the canonical angle is ignored on the grounds of rapid 
phase-mixing in the Langmuir wave, while the particle 
action evolves self-consistently. 

Because these nonlinear, kinetic Langmuir waves arise 
naturally in slowly driven systems, their bulk properties 
can be used to illuminate basic plasma processes and ob- 
tain reduced descriptions of complex phenomena. For 
example, the nonlinear frequency shift of the thermal 
Langmuir resonance is an important quantity in any re- 
duced model of Raman scattering in plasma 0, H, H, [Toll , 
and our results extend those of Morales and O'Neil [11| 
to colder plasmas and larger electrostatic potentials &(z). 
We leave such an implementation in a Langmuir envelope 
code to future work. 

In Sec. |TT] we first present the single particle equations 
of motion and then show that for a weakly driven sys- 
tem, the particles move in an essentially sinusoidal po- 
tential. We proceed by reviewing the relevant pendulum 
dynamics and action-angle coordinates. In Sec. lIIII we in- 
troduce the plasma distribution function, and use a few 
simple, physically motivated assumptions to show that 
a slowly excited plasma is well-represented by a distri- 
bution function that is gaussian in the canonical action 
with a fixed variance. Next, in Sec. lIVl we use Coulomb's 
law and the demands of self-consistency to derive the 
functional relationship between the mean action and the 
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amplitude of the potential. This fully specifies the distri- 
bution function, from which we then extract the natural 
frequency of the BGK-type wave. We compare the mean 
action and frequency of these nonlinear Langmuir waves 
to those obtained from self-consistent particle simulations 
in Sec. [V] for thermal plasmas with 0.4 < k\o < 0.2, 
where Xr> = Vth/^p- Some concluding remarks and pos- 
sible applications are given in Sec. I VII 



II. SINGLE PARTICLE EQUATIONS OF 
MOTION: THE DRIVEN PENDULUM 

In this section we present the single particle equations 
and derive the canonical action-angle variables relevant 
to a weakly-driven plasma wave. These familiar results 
are subsequently used in Sec. IHII to motivate our simpli- 
fied, fluid-like characterization of the electron distribu- 
tion function for a slowly growing Langmuir wave. 

In what follows we ignore transverse variation, as- 
suming that the dominant dynamics are along the lon- 
gitudinal axis z. We ignore the motion of the back- 
ground ions, and furthermore assume that the longitu- 
dinal force on the electrons can be divided into two com- 
ponents: the first is given by an external driving po- 
tential V{z, t) sin(W + kz), which could be, for example, 
a ponderomotive force; the second arises from the self- 
consistent electrostatic potential §(z,t) of the plasma 
electrons. Thus, Newton's equation of motion for the 
longitudinal electron coordinate z(t) is given by 



dt 2 



z(t) 



d_ 



e$(z, t) - V(z, t) sin(wt + kz) 



(2) 



Note that for simplicity we assume the dynamics to be 
nonrelativistic, requiring that the potentials remain suf- 
ficiently small such that |e$| , | | <C m e c 2 ; for a non- 
relativistic external drive V, the electrostatic potential 
will satisfy this relation for all time if the nominal phase 
velocity is much less than the speed of light, (oj/k) 2 <C c 2 . 

To simplify ([2|), we assume that the amplitude of the 
external potential Viz, t) has a slow spatiotemporal vari- 
ation with respect to its carrier phase, so that 



dz 



In|V(js,i)| < k 



d_ 

di 




FIG. 1: Phase space schematic for a monochromatic wave 
overlaid on the results of a self-consistent particle simulation. 
Region I (above the separatrices) consists of the plasma bulk 
making up the wave; region II contains the trapped particles 
between the separatrices; region III contains those particles 
moving too fast to be trapped in the wave. 



Furthermore, since $(z, i) is excited by the slowly vary- 
ing potential V(z, t), the Eikonal conditions © imply 
that the amplitudes 4>n{z, t) and phase shifts £, n {z, t) are 
similarly slowly varying. We assume that the plasma is 
not highly nonlinear, so that the potential is adequately 
described by its first (n = 1) harmonic. We have found 
that this assumption is not overly restrictive: for (f>i as 
large as one-half simulations indicate that 4>2 ~ To^i • 
Finally, we assume that the space-charge force from the 
plasma bulk is dominant, i.e., |V(z,i)| <C |^i(z,t)|. With 
these assumptions, Newton's equation {2} simplifies to 



dt 2 



y-<j>x(z,t) sm[tot + kzj + £i(z, t)]. (5) 



To express (|S|) as a Hamiltonian system appropriate 
for action-angle variables, we introduce the dimensionless 
time r = w p t, the scaled (linear) frequency ujl = w/u> p , 
and the dimensionless coordinates given by the phase in 
the electrostatic wave 9 and its corresponding canonical 
momentum p: 



= cot + kz 



6 



P ; 



- £l = kzj + LU L , 



(6) 



with the over-dot understood to denote the normalized 



time derivative 



dr 



ill 



Defining the frequency shift 



In \V{z, t)\ < to. (3) Slo = ii, © becomes 



The nearly-periodic external drive sets a natural spatial 
length of the driven self-consistent electrostatic potential, 
for which $(z,t) can be Fourier expanded in terms of 
dimensionless Eikonal amplitudes: 



9 =p + 6lo{t) 



P 



(r) s'm9. 



(7) 



The system © can be obtained from the pendulum 
Hamiltonian 



$(z,t) 



"k 2 



Y^ct) n (z,t)cos[n{ujt+kz)+^ n (z,t)}. (4) H(p,6',T 



)=\\p+ 5u(t)] 2 + 20! (r) sin 2 (0/2). (8) 



n=l 



For later convenience, we use a cosine series with the 
Langmuir phase shifts given by £ n (z,t), and the normal- 
ization is chosen such that, for all z and t, 



\<!>n{z,t)\ < 



1 



Here, we include a few important results regarding the 
pendulum Hamiltonian (JSJ. The frozen orbits are de- 
fined as the level sets H(p, 9;t) = H at a fixed time 
r, for which the parameters <pi and Slo are constant and 
the motion is periodic. A representative phase portrait of 
the frozen orbits is in Fig.[TJ superposed on a phase-space 
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snapshot taken from a self-consistent particle simulation. 
Generically, we see that phase space is divided into three 
distinct regions, separated by the trajectories joining the 
hyperbolic fixed points at 9 — ±7r, p = —5lo, for which 
H = 2(f>i. These separatrices separate the rotation mo- 
tion of regions I and III, for which H > 2^>i, from the 
libration about the stable fixed point at 9 — in region 
II, where H < 2</>i. Associated with these frozen orbits, 
there exists a canonical transformation to action-angle 
coordinates Tt(p, 9; r) — > H( J, Vf; t), with the action pro- 
portional to the phase-space area of the frozen orbit: 

J(H-T) = ±-j>d9p(9 1 H-T). 

Using ([5]) , the action is given by the the following integral 
along the frozen orbit: 

J( K; T ) = ^ jd6 ^-(l/K^sin 2 9/2 - Slu(t), (9) 

where we have defined the scaled energy n and the sepa- 
ratrix action X s as 



H 



7T 



(10) 



As is well-known, the integral in can be evaluated in 
terms of complete elliptic integrals. Furthermore, since 
the orbits change their topology at the separatrix, the 
action coordinate is discontinuous there. While no ac- 
tual particle orbits are singular, this discontinuity does 
complicate notation; for this reason, we find it convenient 
to introduce the "frequency-shifted action" X(k) that is 
affinely related to </, and defined according to the phase 
space region in which the trajectory lives. Taking the 
integral in we define 

\k\ > 1 : J(k) = J(k;t) + 5uj(t) =1 s k£(1/k) (11a) 



kl < 1 



I(k) = \[J{k;t)+5u{t)} 



= l s [S{k) + (k 2 - 1) . (lib) 

Here, the complete elliptic integrals of the first and sec- 
ond kind, J4f(n) and (o(k), respectively, are defined in 
the standard way: 



tt/2 



JT(k) 



da 



tt/2 



V 7 ! -k 2 



sin a 



S{k) = da VT - k 2 



sin a. 



Furthermore, the nonlinear period of the pendulum is 
simply calculated using the Hamiltonian relations 

s 2n „ dJ dJ dn 
Tin) = —— = 2ir— = 2vr— — . 

From the definitions (fT0|) and (fTTjl . we have 
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\k\ > 1 : T(k) = —==Jf{l/K) 
\k\ < 1 : T(k) = -%=JT(k). 



(12a) 
(12b) 



Although the transformation (p, 9) — > (J, is canoni- 
cal, we note that the variables (X, ^/) do not form a canon- 
ical pair due to the scaling at the separatrix. Rather 
than use these variables to calculate phase space aver- 
ages, more straightforward analysis is obtained using the 
scaled energy k and the time r. Thus, we conclude by 
relating the coordinates (k, t) to (p,9). Using the defin- 
tions and (jTUJ) , we have 

p + Slu{t) = f- = 2k^[J\- (l/ K 2 )sin 2 (#/2). (13) 
dr v 

Rewriting this expression, we find 
d(9/2) 



k > l 



k < l 



1 - (l/K 2 )sin 2 (0/2) 
da 



K^/fJidr (14a) 
(14b) 



k z sin a 



with sin (0/2) = Ksina. We take the indefinite integral 
of (fl4|) . obtaining 



|k| > 1 : cos(0/2) = cn(l/K, ny/far) 
k| < 1 : cos(6»/2) = dn(/t, ^/far) , 



(15a) 
(15b) 



where, without loss of generality, we have set the origin 
of time to zero, and the functions cn(«, x) and dn(K, x) 
are the Jacobian elliptic functions defined in the usual 
manner via the inverse of the incomplete elliptic integral 
of the first kind: 



for x(k, y) 



dz 



1 



VT 



cos y = crnk, x) = —yn 2 — l + dn 2 (k,x). 

Finally, differentiating (fT5)) and using (fTB")) yields 

|/t| > 1 : p = 2fty^i dn^l/«, Ky^irJ — 6u>(t) (16a) 
|/t| < 1 : p = 2ny/^i cn\K, ^/chrj — 5uj(t). (16b) 



III. THE NONLINEAR DISTRIBUTION 
FUNCTION ACTION ANSATZ 

In this section, we introduce our BGK-type distri- 
bution function ansatz that naturally arises in weakly 
driven plasmas. This results in a dramatic reduction in 
the number of degrees of freedom from the Vlasov equa- 
tion (fTa|) , and therefore can be used as the foundation of 
a simplified, fluid-like model. 

We first motivate our assumed functional form for /, 
showing how simple physical arguments based on parti- 
cle phase mixing in the slowly changing electrostatic field 
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governed by Vlasov evolution lead to a very simple dis- 
tribution: gaussian in the canonical action J with slowly 
evolving mean and fixed variance. We include an exam- 
ple from a slowly driven, self-consistent particle simula- 
tion that naturally evolves in this way, and compare it 
to the more familiar velocity-space distributions. Finally, 
we give the distribution function in an explicit form from 
which phase-space averages can be computed. 



A. Phase-mixing in the slowly-growing wave 

The central assumption for our distribution function 
action ansatz is that the Langmuir wave amplitude and 
frequency are slowly evolving, meaning that the parame- 
ters 4>i{t) and Slo(t) of the pendulum Hamiltonian ([8]) do 
not vary appreciably over the period T . For the nearly all 
electrons, i.e., all except the exponentially few in a nar- 
row range about the separatrix, the condition of "slow- 
ness" can be written as 



1 d 



'dr 



line 



1 



hi dr 



|ln<M 



e< 1. 



(17) 



As previously noted, while these conditions are rem- 
iniscent of adiabatic motion, we do not explicitly in- 
voke adiabaticity; rather, we use the slowness condition 
(I17p to justify our assumption that the distribution func- 
tion remains essentially uniform in the canonical angle 
\1/ throughout its evolution. Far from the separatrices, it 
is clear that the slowness conditions (fT7|) imply that the 
electrons make many oscillations before the parameters 
of the wave significantly change, so that a set of these 
particles that is initially uniform in canonical angle re- 
mains so under evolution by ([5]). As particles approach 
the separatrix where k — > 1 and the nonlinear period 
diverges ~ In 1 1 — k 2 | , such a simple argument breaks 
down. In this case, we invoke the results of phase evolu- 
tion by Cary and Skodje [HI, EH and Elskens and Escande 
[l4| , obtained by analyzing the near-separatrix motion in 
slowly evolving systems. 

For most of the particles, the canonical angle is 
mapped smoothly through the separatrix. In a naive pic- 
ture, the strip of particles in region I with the same action 
J and spread over < ^ < 2n is mapped across the sep- 
aratrix to the strip from < '3/ < ir that then rotates 
in region II. As shown in [l2l . Il3l . [Tij , this picture is es- 
sentially correct to O(e), excluding the exponentially few 
0(e~ 1 / e /e) particles that pass very close to the hyper- 
bolic fixed point. Since these particles can spend an ar- 
bitrarily long time tracing the stable manifold, they lead 
to long, diffuse phase-space tendrils. Neglecting these 
particles, to each action is associated a strip of parti- 
cles that is mapped to one-half the canonical angle upon 
crossing the separatrix. This proceeds with each succes- 
sive action strip, with each one displaced from the next 
by a relative phase in the canonical angle. The relative 
phase between increasing action strips increases up to 
2ir, for which the action has increased by 0(e) (see 14 1 



for a detailed discussion). In this way, we argue that 
the distribution remains phase mixed to 0(e) even when 
crossing the separatrix, provided only that the slowness 
conditions (TTTj) are met. 

Assuming that the electrostatic wave amplitude and 
phase velocity are slowly varying, the previous arguments 
indicate that the distribution function remains uniform 
in the canonical angle throughout the range < 'J < 2tt, 
resulting in the following simplification: 



/(■/,*;r) -£/(./; r). 



(18) 



To make further progress, we also assume that the dis- 
tribution /(J;r) is well-represented by its first two mo- 
ments in canonical action. This is trivially true if </>i = 0, 
for which v oc J and a Maxwellian plasma is gaussian in 
action. In the general case, this assumption is similar in 
spirit to a warm-fluid theory (see, e.g., [HI, EH]), for which 
asymptotic solutions are obtained by systematically ne- 
glecting the (presumably small) third- and higher-order 
moments in v. By using moments in J, however, our 
model will self-consistently include the effects of trapped 
particles on the Langmuir wave, as do the models of Hol- 
loway, Doming and Buchanan [y, [l7| ; on the other hand, 
unlike the theories of [HI, [III, ours is non-relativistic. 
Retaining only the first two moments in the canonical 
action is equivalent to prescribing / to be gaussian in J: 



f(J;r) 



1 



a(T)V2Tr 



exp 



[J-J(r)f 
2<j(t) 2 



(19) 



We can simplify (TIT))) in the following manner. The 
Vlasov equation (|la[) permits an infinite number of con- 
servation laws (Casimirs), of which the entropy is one of 
particular physical significance: 



d 

~d7 



fdJd* /(J,*;r)ln/(J,*;r) 



0. 



(20) 



Using the gaussian-in- action distribution (|19p . the con- 
servation of entropy (f20|) implies that 



dr 



In cr(r) = 0, 



(21) 



i.e., that the gaussian width is fixed throughout the evo- 
lution. Some care must be taken to apply this result 
across the separatrix where the action J is not contin- 
uous; note that this is an issue of coordinates, not dy- 
namics. Since the coordinate J doubles across the sepa- 
ratrix, the width in J similarly doubles. To be more pre- 
cise, the phase space region {[Ji, J\ + SJ], [0,27r)} just 
above (or below) the separatrix corresponds to the equal 
area region between the separatrices over the intervals 
{[2Ji,2(Ji +SJ)], [0,tt)}. Thus, in terms of the scaled 
actions I(n) defined in (fTTj) . the width a is the same in 
all regions. Using the results ([!§)) and (|2"Tj) , we arrive at 
the following form for the distribution function 



/(«;t) 



i l 

2tt cfJ2tt 



exp{-^[X( K )-J(r)] 2 }. (22) 
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FIG. 2: Wavelength averaged distribution functions for three 
values of the electrostatic potential (pi using feAn = 0.3. 
In velocity-space, (a) shows the flattening of f(v) about the 
phase velocity kv/u> p ~ 1.16. Meanwhile, (b) demonstrates 
that the distribution in canonical action remains nearly gaus- 
sian throughout the excitation of the plasma wave. 



Note that as defined the initial width is a measure of the 
temperature such that a = k\p. 



To illustrate the distributions associated with ([22)1 , we 
have performed a number of single-wavelength particle 
simulations, described in greater detail in Sec. [V] that 
solve the periodic Vlasov-Poisson system. We include 
representative results in Fig. [2], obtained with a drive po- 
tential V = 0.01 and initial a = 0.3. In Fig. HJa) we see 
the characteristic flattening of f(v) near the phase veloc- 
ity that is associated with particle trapping. For the same 
values of <f>x , Fig. [2^b) demonstrates that the distribution 
in J (integrated over 'F) remains nearly gaussian in the 
canonical action. Furthermore, except for the slight os- 
cillations near the resonance region, /(J) has a constant 
variance a 2 and a slightly increasing mean (from (Jl). 
Note that Fig. [2] uses the numerically calculated J from 
the exact simulation potential, which we find only differs 
from the pendulum results very near the separatrix. 



B. Phase space averages 

Here, we combine the distribution function (|22[) with 
our assumptions on the slowly- varying nature of X to ar- 
rive at a convenient expression of phase space averages 
in terms of the variables k and r. First, we note that for 
fixed parameters </>i, 5u>, there exists a canonical trans- 
formation of Ti : (p, 9) — ► (H, t) for which the evolution 
parameter is the coordinate 9. Since the transformation 
is canonical, the Jacobian is unity and we have the fol- 
lowing relation between the integration measures: 

dH 

dpd9 = dH dr = ——dndr = 4d>in dndr. (23) 
dn 

Thus, phase space averages can be written as 



(X{p,e))= dp d9 f(p,6;T)X(p, 



oo T 

= J dn 40i « J dr f(n;T)X(n, r). (24) 

-oo 

Assuming that the average action does not change signifi- 
cantly during one period, i.e., that the conditions (|17[) are 
met, we can take the nearly constant distribution /(k; t) 
through the integral over r. Hence, in subsequent calcu- 
lations we suppress the dependence of / on r. To further 
simplify notation and integration limits, we assign the 
following expressions for the distribution in k appropri- 
ate to the various phase space regions: 



/iu(k) 
/n(«) 



^^exp{-^[l(.)-l] 2 } (25a) 
^ (1/K) exp{^[l(.) + I] 2 } (25b) 



(25c) 



oV27T 



,KJZ (K r j f -,21 

(Tv Zn *■ J 



The definitions (f2"5)) arise by absorbing the factor of 40ik 
into the / of (f2"2"| . appropriately changing the sign of k, 
and multiplying by the period T from (fT2|) ; these con- 
ventions imply that 



dn [£(«) + /m(«)] + / dn [/„(«) + /+(«)] = 1. 



Using the definitions (|25|) in ([24|) . and remembering the 
chosen sign conventions for k and the scaling by T, phase 
space averages are computed with the integral expression 



(X) = Idnldr + fm{K)Z%% 





1 T 

dn J dr 

o 



(26) 
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IV. NONLINEAR SELF-CONSISTENCY: 
PARAMETERIZING THE DISTRIBUTION 

Now that we have presented the distribution function 
ansatz and developed the necessary pendulum machin- 
ery, we turn to calculating expressions for the mean ac- 
tion X and the frequency shift Su>. We do this by impos- 
ing the constraints of self-consistency: the assumed dis- 
tribution function must also satisfy Maxwell's equations. 
After obtaining X and 8u> in terms of integrals over the 
modulus of the elliptic functions, we will derive small- 
amplitude, analytic expressions suitable to comparisons 
with previous published results. In the next section we 
compare these theoretical results to numerical examples 
from a self-consistent particle code. 



A. Self-consistency and Maxwell's equations 

First, our nonlinear distribution function must give rise 
to a charge separation commensurate with the electro- 
static potential fa[z,t). This is summarized by the Pois- 
son equation: 



d 2 $ 
Ik? 



Aire / dv f(v, z; t) — 47reni(z). 



(27) 



Fourier transforming (|27ll , assuming the ions are station- 
ary [rii{z,t) = no], and using the dimensionless Fourier 
expansion of the potential ([4J, the Poisson equation for 
the lowest harmonic of the potential can be written as 



= 



1 + — (cos 9(k,t)) 
fa 



fa = e(X;fa)fa, (28) 



where we have introduced the "nonlinear dielectric" 
e(X;fay From (|2"5|) , we see that for our assumed distri- 
bution to support a nontrivial potential fa , the nonlinear 
dielectric must vanish. This is an implicit relation for the 
mean action I in terms of the potential amplitude, i.e., 
given a potential fa, the requirement e(X; fa*} = estab- 
lishes the appropriate value of X. To determine the non- 
linear dielectric, we use our assumed distribution func- 
tion to calculate the phase space average of cos 9(k, t) as 
indicated in (|28|) . We define the variables x = K\ffar, 
y = \ffar, and use the pendulum identity (TT5)) and the 
phase space averaging (f2l>|) , to arrive at 



JT(1/k) 



(costf) = dn [/[(«) + /n(«)] / dx 



2cn 2 (l/ K ,a;) - 1 



JT(/s) 



+ dK [fu(*) + /£(*)] / dy 



JT(1/k) 



2dn 2 ( K ,y) - 1 



The integrals over x and y can be taken analytically; 
using the integral tables of Gradshteyn and Ryzhik [lj|, 



the nonlinear dielectric is given by 

oo 

e(X,fa)=l + -^-JdK + /iii(k)] 

1 

'2k 2 S{1/k) 



JT(1//s) 



+ - dn [/+(«) + fn(*)] 



1 - 2k 2 



(29) 



2£{k) 



- 1 



Requiring 129|) to vanish gives an implicit relationship 
between fa and I; to obtain an explicit expression for the 
mean frequency-shifted action in terms of the potential, 
we expand the nonlinear dielectric for small changes in X. 
For an initially Maxwellian plasma with zero electrostatic 
field, the mean action in the moving frame is equal to the 
linear frequency ujl, and we have the expansion 



d 

= e(u)L, <t>i) + wl) Tjj e(X, fa 



(30) 



Rewriting the Taylor expansion (l30|) . the change in the 
mean action in terms of potential amplitude is 



X - lo l w - 



e(uL,fa) 



(31) 



Equation (|3ip determines the mean frequency-shifted ac- 
tion X required to support a potential of amplitude fa . In 
previous works (e.g., [1, HH), expressions similar to ((3"T|) 
were interpreted as the frequency shift of the wave. We 
will see that in the small fa limit, the change in X equals 
the change in the frequency. This is because the parti- 
cle action J is essentially constant in this case, so that 
X oc p implies that a change in X is due to a decrease 
in the phase velocity of the wave at fixed k [i.e., given 
by Sll>(t)]. As the wave amplitude becomes appreciable, 
however, the plasma bulk becomes excited and the in- 
crease in canonical action J can no longer be neglected, 
so that 5lu{t) < X{t) —ujl and ([3Tj) is numerically larger 
than the (generically negative) frequency shift. 

To complete the characterization of the action distri- 
bution function, we calculate an expression for the fre- 
quency shift as a function of fa and the mean action X. 
To determine Sco, we use the fact that the plasma tends 
to set up a return current to erase any long-range elec- 
tric fields. This is related to the fact that in 1-D with 
immobile ions, the electrons gain no net momentum. To 
make this more explicit, we consider the 1-D Coulomb 
gauge condition: V • A = -§^(z ■ A) =0, implying that 
the vector potential can be taken to be purely transverse, 
i.e., A ■ z = 0. In this case, the longitudinal component 
of the Ampere-Maxwell equation is given by 



i d 2 



dtdz 



$(z, t) = 47re / dv f(v, z; t) v 



(32) 
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Integrating over one period in z, we have 

OO 7V 

[*(*,*)_ 3 (_*,*)] =4ne jdv Jdzfv. (33) 



— OO — 7T 



Since we assume the potential $(z. t) to be slowly vary- 
ing, the expression ()33|1 approximately vanishes. Note 
that this is exact in the limit of a time- independent, non- 
linear mode (similar to BGK), and expresses the fact that 
the plasma electrons carry no net momentum [l9l | . Using 
the fact that the velocity v cx (p — u>l), we have 



(p(k,t)) -uj l =0. 



(34) 



To evaluate this phase space average, we use the pendu- 
lum formula (|16[) . which gives the momentum in terms 
of Jacobi elliptic functions and the frequency shift 5u> (r). 
Since the averaged velocity of the trapped particles in the 
moving frame is zero [the integral of cn(«, x) from to 
4J?T(k) vanishes], we have 



ujl + Soj(t) = I dn 



Mk) - /iii (k) 
JT(1/k) 



JT(1/k) 



J dx dn(l/K, x) 



We take the integral over x = n^fcfiT using Gradshteyn 
and Ryzhik [lH, so that 



5u> (r) — dK 



3lK 



e 2o- : 



' r i a K^(l/ K )-Xl 2 



#1^ i [i sK(? (i/ K)+ i] 2 _ 



(35) 
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B. The mean action 1 and frequency shift 5u>: 
linear and small amplitude limits 

In this section we study the linear and small amplitude 
limits of X and Slo, and compare them to previous results. 
First, we present the linear limit of (j3"Tj) . for which <\>\ — » 
0. In this limit, the mean action is that corresponding 
to the phase velocity of the infinitesimal wave, so that 
X — * u>l. Since the denominator in (|31[) is well-behaved, 
the linear limit is characterized by 



lim e(w L ,0i) = 0. 

<Ai — >0 



(36) 



For clarity, we reserve the cumbersome calculations used 
to evaluate (|36p for the Appendix. In the Appendix, we 
show that our assumed distribution gives a concrete pre- 
scription for the usual pole occurring when the particle 
velocity matches that of the wave phase velocity. Denot- 
ing the principal value by 3 s , from (|A.3|) we have 



lim £(wl,0i) 
01— >o 



1 



LU L /<7 



/2vr 



9> I dx 



uj L /a 



Setting the linear dielectric (|37| to zero yields the plasma 
dispersion relation as found by Vlasov [20| , resulting in a 
purely real natural frequency. Physically, this lack of lin- 
ear Landau damping arises because we have assumed that 
the distribution is completely phase mixed. As shown 
by O'Neil such phase mixing causes linear Landau 
damping to be a transient effect that itself damps away 
on the bounce time scale ~ 1/(w p a/0i")- Although the 
bounce time diverges as <j)\ — > 0, we maintain that our 
analysis and the dispersion relation (|37[) applies to fi- 
nite amplitude waves after several bounce periods have 
passed. In this case, the distribution is nearly uniform in 
canonical angle and Landau damping has been "washed 
out" . The nonlinear fate of such Langmuir waves is gen- 
erally considered to be a BGK wave; as discussed in 
and [l7| . the dispersion relation of small amplitude, si- 
nusoidal BGK waves is that of Vlasov, and is identical to 
([57)1 derived here. 

For 4>\ small but not infinitesimal, we Taylor expand 
the dielectric ([29]) in the Appendix, yielding 



1.089a 
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(38) 



A similar expression for the small amplitude dielectric 
has been derived by a number of authors, although there 
is some variation in the 0(1) coefficient replacing our 
1.089. Interestingly, our coefficient is precisely that of 
Dewar (22j , who calculated the frequency shift assuming 
a small but finite sinusoidal wave that is adiabatically 
excited; other calculations in a similar regime have ob- 
tained values of 1.41 (Manheimer and Flynn [23|), 1.76 
(Rose and Russell @) and 1.60 (Barnes [24j]). It should 
be noted that these differ slightly from the coefficient of 
1.63 calculated by Morales and O'Neil [ll[ and separately 
by Dewar [22l | for an instantaneously excited wave (i.e., 
the initial value problem). The majority of these authors 
then use this to determine the nonlinear frequency shift 
using an expression similar to (|3 1 [) : we will see that this 
yields reasonable results provided that <f>\ < a 2 . 

To complete the small amplitude analysis of the change 
in the average action ([3"Tj) , we differentiate the linear dis- 
persion relation (|37p . obtaining -t^-e(ujl, <f>\) in the small 
amplitude limit. Using ([3"Tjl and ([3"8]) , we find that the 
change in the mean action is expressed for small values 
of 0i by 



T(t) - u L 



-1.089 wi,-. 



-1.089a 



CM 



2 In 2 



2ttct 



1 



-/"(J) J=0 , (39) 

J = UJ L 



with the distribution f(J) given by (Q3 

To conclude this section, we give the difference between 
the frequency shift 5lo and the shift in action X — u>l,, 
showing that these differ as the plasma wave amplitude 
(37) <f>\ grows. In the Appendix, we see that this difference is 
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FIG. 3: Frequency shift for three different temperatures: fcAr> = 0.4 (a), /cAd = 0.3 (b), and fcAr> = 0.2 (c). The points are from 
simulation results, with error bars indicating two standard deviations in the measured data. We see that the theoretical value 
of 5ui from the action distribution (solid line) agrees quite well with the simulation results, which are only closely represented 
by the Dewar/O'Neil theory (|39|1 for small values of 4>i- To calculate our theoretical frequency shift, we evaluate the numerical 
integral (135[) with the mean X evaluated numerically using the approximate expression (131[) . 



given by 

duj w [1 - UJ L ) ~ T 7=^—0! + 

V. COMPARISON TO SELF-CONSISTENT 
PARTICLE SIMULATIONS 

In this section, we compare our theoretical results for 
the properties of the slowly-driven nonlinear Langmuir 
waves with those obtained from particle simulations. Be- 
fore discussing these examples, we make a few comments 
on the numerical methods. In these single-wavelength 
simulations, we numerically solve the equations of motion 
for the electrons and the electric field, driven by an exter- 
nal force. For a wavelength with TV electrons, the electron 
with coordinate Q = k^Zj experiences the combined self- 
consistent electrostatic force and prescribed drive, giving 
rise to the following equation of motion: 

d 2 M 1 N 2 

^(r) = £ n S m sin K,(T) - mfc(r)] 

m=l 1=1 v > 

+ V(t) cos(uj l t + Cj) - £o(t), 

where we have expanded the electrostatic potential in M 
harmonics, each of which is a sum over the N macro- 
particles. This is a standard technique of the free elec- 
tron laser community [2 51, although here we have also 
retained the DC field £q [26|], to be calculated using the 
longitudinal component of the Ampere-Maxwell law: 

, 1 N , 

For the examples shown here, we use N sa 10 6 particles 
and M « 32 harmonics. We solve the system (|40M41[) for 
a given drive potential V(r) using a symplectic operator- 
splitting scheme that is second-order accurate. 



To compare the simulation results to our theory, we 
slowly turn on the ponderomotive drive, ramping the 
electrostatic field to a chosen amplitude, at which time we 
slowly turn the drive off. By taking the Hilbert transform 
of the potential we obtain the total frequency (l>l + 8u>), 
from which we extract the frequency shift for a given am- 
plitude <f>\. These results are shown in Fig. [31 where we 
compare the frequency shift extracted from simulation 
to the well-known result of Morales and O'Neil [ll[ and 
Dewar [22] , Sui ~ vW/" (X) > an d to our theory for three 
different values of fcA^, = a: 0.4 (a), 0.3 (b), and 0.2 (c). 
For the Dewar/O'Neil theory, we use equation (|39[) which 
is of the form first derived by Morales and O'Neil [TU ] . 
but with the O(l) pre-factor of Dewar [221 ]. We evaluate 
the action ansatz value of 5u by numerically integrating 
with the average action J given by the relations 
(US]) and d3U). 

As we can see in Fig. [31 the Dewar/O'Neil theory 
agrees with our results for small values of <p±, but de- 
viate at larger values of the potential. Furthermore, 
our nonlinear theory captures both the qualitative and 
quantitative features of Su> seen in the particle simula- 
tions over the entire range <pi . The discrepancy between 
theories becomes particularly clear for colder plasmas; 
for k\o = 0.2, (c) shows that the Dewar/O'Neil the- 
ory is only applicable when there is negligible frequency 
shift, while the action ansatz yields reasonable qualita- 
tive agreement for all amplitudes. We do note, however, 
that the quantitative agreement is worse than that ob- 
served for warmer plasmas; we speculate that this is due 
to the significant harmonic content of these highly non- 
linear Langmuir waves. 

The range of <j>i over which 5u was measured in Fig. [3] 
includes all electrostatic amplitudes that were attained 
via resonant excitation with the drive amplitude V — 
0.01 and frequency equal to wl- Further driving of the 
plasma results in a ringing of <f>\ that we interpret as 
resulting from the de-tuning of the nonlinear Langmuir 
wave from the external drive. 
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0.1 0.2 0.3 



FIG. 4: Change in the canonical action J for three different 
temperatures: fcAo = 0.4, fcAzj = 0.3, and fcAu = 0.2. The 
points are obtained from the simulation while the curves are 
numerical evaluations of J using J = T — Suj. As in Fig. [3] the 
range of <f>\ includes all amplitudes attained with a V = 0.01 
resonant drive; we interpret the saturation of <j>\ to be caused 
by de-tuning from the drive. 



VI. CONCLUSION 

We have presented a nonlinear electron distribution 
function that naturally arises for slowly driven Lang- 
muir waves, and is parameterized by the amplitude of the 
electrostatic potential. The distribution function is de- 
scribed by a gaussian in the canonical action whose mean 
varies self-consistently according to the slowly evolving 
potential, and whose variance remains fixed at the initial 
plasma temperature. We then determined the nonlinear 
frequency shift of the Langmuir wave (|35[) , which agrees 
well with full particle simulations. This frequency shift 
could be used in a reduced, fluid-like model for driven 
plasmas, while the asymptotic distribution function (|22|) 
hints at a simplified model of nonlinear Landau damping 
via the dynamical process of phase-mixing. Present re- 
search aims to use these results as the kinetic foundation 
of an extended three-wave model of Raman backscatter 
in a thermal plasma. 

Acknowledgments 



Finally, we conclude this section with simulation re- 
sults regarding the measured change in the mean canon- 
ical action J. To measure this, we use the extracted 
frequency shift 5u> to determine the energy of each par- 
ticle using the total potential (including harmonics): 
Hj = ^(pj+Sui) 2 — (f)((j). Solving this for the momentum, 
we numerically integrate the exact frequency-shifted ac- 
tion for each particle, and obtain the average J via 



J = 



1 E- 



dC J 2 [Hj + 0(C)] - Slu. 



(42) 



In the usual adiabatic approximations, the particle action 
is conserved so that (|4"2"|) is constant. However, our theory 
predicts J to increase as the potential grows. We inter- 
pret this as resulting from the manner in which we excite 
the plasma: since the plasma is resonantly driven, the 
time-scale separation between the drive and the plasma 
response can be vanishingly small even for slowly chang- 
ing amplitudes. Thus, the particle action need not be an 
adiabatic invariant. 

To further validate this claim, we present the simula- 
tion results for the change in the average canonical action 
J in Fig. 21 and compare this to our theory for the tem- 
peratures of Fig. [3] We see remarkable agreement for 
temperatures such that feAn = 0.4, 0.3, 0.2. First, this 
indicates that the canonical action of the particles is not 
conserved in these slowly-driven problems. Rather, we 
find a well-described relation between the electrostatic 
potential and the change in the average action. In the 
small amplitude limit, the change in J scales as <f>\^ 2 as 
shown in the Appendix (|A.8[) . while numerical evidence 
indicates that it scales roughly as <j)\ for <fii > a 2 . 
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APPENDIX: LINEAR AND SMALL AMPLITUDE 
INTEGRALS FOR I, Slj 

To evaluate the linear and small amplitude limits of 
the dielectric e{ujL,(j>i), we first integrate lp2"9")) by parts. 
The boundary terms at k — 0, ±oo vanish, while those 
at k — 1 cancel, resulting in the following formula for the 
nonlinear dielectric: 



e(wi,</>i) 



dn 



h(n) 
~T7 



(Ik 



where 

h{n) 



32 



37T 2 cr 2 

32 

3lT 2 <J 2 



+/m(/c) T{k) + lo l | 

{/n( K ) J ( K ) ~ u l 
+/+(«) [Zto+Wi]}, 

(2k 3 - k)S{1/k) - 2(k 3 - k)JT(1/k) 
(2k 2 - 1)<?(k) - (k 2 - l)jr(re) 



(A.l) 
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First, we calculate (IA. 1|) in the linear limit, for which 
4>i — ► 0. In this case, the integrals over the trapped 
electrons (for which < k < 1) make no contribution, 
while the remaining integrand vanishes exponentially for 
small K. Thus, we consistently take n 3> 1, for which we 
have 



h{n) 



- + 0(1/k 3 



e 2a 



2ny 



X(k) =T s K§(\jK) 

Defining the variables 

ax = 2nyQ)i + u>l ay 



27TCT 3 

0(1). 



(A.2a) 
(A.2b) 



and using the large k relations (|A.2|) . we find that the 
linear relation lim = is given by 



= 1 + — 

a 1 



lim 



dx ■ 



-x 2 /2 



2na 2 S—>o | J x-uj L /a 



(A.3) 



dy 



-y 2 /2 



y - uj L /a \ ° 



where 8 = 2v^i- As we can see, the assumed distribu- 
tion function yields a prescription for treating the pole 
at x = uj^/a: the symmetric limit is merely the prin- 
cipal value. Note that this is the standard pole occur- 
ring when the particle velocity equals the phase velocity 
of the wave (i.e., when the particle action equals that of 
the separatrix defined by the infinitesimal wave) , and the 
symmetric limit results in the Vlasov-type dispersion re- 
lation (|37p . Furthermore, differentiation yields the quan- 
tity 0) relevant for the calculation of the change 
in the mean action. Using the dispersion relation (|37[) . 
we find 



lim — — 

01^0 0LU L 



e{uj L ,(j>i 



1 



LU L a z 



M. -I- a 2 ) 



(A.4) 



We continue by calculating the change in X (|3T|) in- 
duced by the near-resonant particles, assuming the am- 
plitude of the potential <j>i is small. To make the inte- 
grals defining e(u>l,<Pi) manifestly convergent, we start 
by first "rewriting the 1" in the expression for the dielec- 
tric (jA.lj) . Assuming the linear dispersion relation ((37|) 
is satisfied, in terms of the energy k we have 



\ dn 



2ka 



e 2 
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(A.5) 



By appropriate choice of signs for k, we can express (|A.5[) 
as a sum of integrals whose limits are such that 1 < k < 



oo or < k < 1, that we then use to replace the 1 in the 
nonlinear dielectric (|A.1|) . Thus, the nonlinear dielectric 
is given by 



e(lu l ,4>i) = I < Ik 
1 



2(k) - uj l h(K)fi(n) 
dn X{n) +lul h(n) fin(K) 
1(k) -cj l g(«)/n(«) 



(A.6) 



dn 



1(k)+lu l q(K)f+(K) + l, 



where the number 1 is to be interpreted as the sum of in- 
tegrals from (|A.5[) . In this manner, the expression (|A.6|) 
is perfectly well-defined in the small amplitude limit, and 
this limit is simple to calculate numerically. Taylor ex- 
panding the integrals with 1 < k < oo, we have 



f>i - — - I di; 



3tt 2 



37T a 

— -h{K)Kg{l/R)Jf{l/K) 



aV2n 



-1.50a 



D -^l/2a 2 
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while the last two integrals of (|A.6|) yield 



64k 



l^x I dK \ 1 - g-g g(«) JT(k) [ S{k) + {k 2 - 1) JT(k)] 



x 4 



-u> 2 L /2a 2 



2.59a 



\/2° 2 



2na 3 



Adding these contributions, we find that the small am- 
plitude behavior of the nonlinear dielectric is 



lim e{u L ,<t>i) = 1.O89V0 

01^0 



1 



u> 2 L /2a 2 



2na 3 



(A.7) 



Finally, we calculate the difference between the fre- 
quency shifted average action and the frequency, namely, 
the average action J. Using the definition X = (X) and 
the expression for the frequency shift ([35)) , we find that 



[u>L + SuA = dn 



X(k) 



Jf(l//s) 



+ JdKX(n) («)-/£(«)]. 
n 

It can be shown that both integrals vanish in the limit 
<f>i — > 0. Taylor expanding the integrals for small values 
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of the potential <j)\ , we find that 

32 u L e-»y 2 °\ V2 




—5- <?(!/«) JT(1/k) - k z 



+ Jdn^J(f{K) [<?(«) + (k 2 - 1) JT(k)] 

(A.; 



64 w L e^/ 2 * 2 



971^ aV2^ 1 



— 3/2 

Thus, we see that the average action J grows as e\ f° r 
0i small. For Langmuir wave amplitudes such that these 
terms can be neglected, the frequency shift is equal to the 
change in the frequency shifted action, and the canonical 
action is approximately conserved. 
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